00100		SUBROUTINE LPC1(FFY,CF,RO,SPT,NPTS,M,NSP)
00200	C	THE PARAMETERS IN THE LIST ARE DEFINED AS FOLLOWS:
00300	C	IFFY←THE INPUT DATA ARRAY OF INTEGERS...A TIME SERIES
00400	C	NPTS←THE NUMBER OF POINTS IN THE TIME SERIES
00500	C	CF←  THE ARRAY IN WHICH THE FILTER COEFFICIENTS
00600	C		ARE TO BE RETURNED.
00700	C	M←   THE NUMBER OF FILTER COEFFICIENTS DESIRED.
00800	C	RO←  COMPUTED FILTER SCALE FACTOR.
00900	C	ERRN←NORMALIZED ERROR VALUE.
01000	C	ERR← ARRAY IN WHICH THE SERIES OF ERROR IS RETURNED.
01100	C	SPT← THE ARRAY IN WHICH THE REAL SPECTRUM IS RETURNED.
01200	C	NSP← THE NUMBER OF POINTS IN THE REAL SPECTRUM.
01300	C	ISSW← AN INTEGER SWITCH USED TO SELECT THE FOLLOWING
01400	C		MODES OF OPERATION: **** NOW REDUNDENT ****
01500	C		1: CALCULATE AND RETURN ALL VALUES IN LIST.
01600	C		2: CALCULATE RO, ERRN, AND CF.
01700	C		3: CALCULATE RO, ERRN, CF, AND ERR.
01800	C		4: CALCULATE RO, ERRN, CF, AND SPT.
01900	C		5: USING CF AND M AS INPUT CALCULATE ERR AND SPT.
02000	C		6:   "    "  "  "  "   "       "     ERR.
02100	C		7:   "    "  "  "  "   "       "     SPT.
02200	C		1X: FOR ANY X ABOVE IMPLYING IFFY INPUT,
02300	C		    SUPPRESS HAMMING WINDOW APPLICATION.
02400		DIMENSION FFY(1),CF(1),ERR(1),SPT(1),X(512),Y(512)
02500		TPI=2.*355./113.
04100	10	NP=NPTS-1
04200	C	CALCULATE THE DIFFERENCE SIGNAL
04300		DO 11 J=1,NP
04400	 11	X(J)=FFY(J)
05300	C	CALCULATE AUTOCORRELATION COEFFICIENTS AND SOLVE
05400	C	FOR INVERSE FILTER COEFFICIENTS
05410	100	DO 105 J=1,NPTS
05420	105	Y(J)=X(J)
05500		CALL INVFL(Y,NPTS,M,RO,ERRN)
05600	C	Y ARRAY RETURNS 1.,A1,A2,...,AM,0.,....0.
05700		MM=M+1
05800		DO 125 J=1,MM
05900	125	CF(J)=Y(J)
07200	300	DO 303 J=1,512
07300		X(J)=0.
07400	303	Y(J)=0.
07500		MM=M+1
07600		DO 310 J=1,MM
07700	310	X(J)=CF(J)
07800	C	SHUFFLE DATA FOR CALC OF REAL TRANSFORM
07900		MM2=MM/2 +1
08000			DO 320 K=1,NSP
08100			K1=2*K-1
08200			K2=2*K
08300			X(K)=X(K1)
08400	320		Y(K)=X(K2)
08500	 	NMOD=1
08600	  323	IF (2**NMOD.GE.NSP) GO TO 326
08700		NMOD=NMOD+1
08800		GO TO 323
08900	 326	IF ((2*MM2).LT.MM) MM2=MM2+1
09000		MOD=1
09100	 330	IF (2**MOD.GE.MM2) GO TO 340
09200		MOD=MOD+1
09300		GO TO 330
09310	C	NMOD=LOG(NSP)
09320	C	IF (2**NMOD.LT.NSP)  NMOD=NMOD+1
09330	C	MOD = LOG(MM2)
09340	C	IF (2**MOD.LT.NSP)  MOD=MOD+1
09400	340	CALL FFTP(X,Y,NMOD,MOD)
09500		CALL RBITS(X,Y,NMOD)
09600		CALL USCRM(X,Y,NMOD)
09700	C	CALCULATE LOG OF RECIPROCAL INVERSE FILTER SPECTRUM
09800			DO 350 J=1,NSP
09900	350		SPT(J)=(X(J)*X(J)+Y(J)*Y(J))
10000	C		DO 360 J=1,NSP
10100	C 360		SPT(J)=-10.*ALOG10(X(J)*X(J)+Y(J)*Y(J)+1.)
10200		RETURN
10300		END
     

00100		SUBROUTINE INVFL(X,NP,M,RO,ERRN)
00200		DIMENSION X(512)
00300		DIMENSION F(50),TF(50),A(50),R(50)
00400	C	CALCULATION OF M+1 LENGTH AUTOCORRELATION SEQUENCE
00500	C	FROM THE INVERSE FILTER INPUT SEQUENCE
00600		MP1=M+1
00700			DO 11 JJ=1,MP1
00800			J=JJ-1
00900			MMJ=NP-J
01000			SS=0.
01100				DO 10 I=1,MMJ
01200				IPJ=I+J
01300	10			SS=SS+X(I)*X(IPJ)
01400	11		R(JJ)=SS
01500		RO=R(1)
01600	C	FORTRAN IMPLEMENTATION OF LEVINSON'S METHOD FOR
01700	C	SOLVING THE AUTOCORRELATION MATRIX
01800		F(1)=1.
01900		ALPHA=R(1)
02000		BETA=R(2)
02100		A(1)=-R(2)/R(1)
02200		GAMMA=A(1)*R(2)
02300			DO 1 N=2,M
02400			NM1=N-1
02500			C=-BETA/ALPHA
02600			IF (N-2) 2,2,3
02700	3			DO 4 J=2,NM1
02800				NN=N-J+1
02900	4			TF(J)=F(J)+C*F(NN)
03000				DO 5 J=2,NM1
03100	5			F(J)=TF(J)
03200	2		F(N)=C*F(1)
03300			ALPHA=ALPHA+C*BETA
03400			BETA=0.
03500				DO 6 J=1,N
03600				NN=N-J+2
03700	6			BETA=BETA+F(J)*R(NN)
03800			Q=-(R(N+1)+GAMMA)/ALPHA
03900				DO 7 J=1,NM1
04000				NN=N-J+1
04100	7			A(J)=A(J)+Q*F(NN)
04200			A(N)=Q*F(1)
04300			GAMMA=0.
04400				DO 8 J=1,N
04500				NN=N-J+2
04600	8			GAMMA=GAMMA+A(J)*R(NN)
04700	1		CONTINUE
04800	C	CALCULATE NORMALIZED ERROR ERRN
04900		SM=0.
05000			DO 77 J=1,M
05100	77		SM=SM+A(J)*R(J+1)
05200		ERRN=1.+SM/R(1)
05300	C	PLACE COEFFICIENTS IN PROPER LOCATIONS OF THE
05400	C	REAL VECTOR X FOR FFT-ING TO GET INVERSE SPECTRUM
05500			DO 51 J=1,M
05600	51		X(J+1)=A(J)
05700		X(1)=1.
05800		MP1=MP1+1
05900			DO 123 J=MP1,512
06000	123		X(J)=0.
06100		RETURN
06200		END
     

00100		SUBROUTINE FFTP(X,Y,M,L)
00200		DIMENSION X(512),Y(512)
00300		N=2**M
00400		L2=2**L
00500		        DO 1 LO=1,M
00600		LMX=2**(M-LO)
00700		LMM=LMX
00800		LIX=2*LMX
00900		SCL=6.283185/LIX
01000	C     TEST FOR PRUNING
01100		IF (LO-M+L) 20,30,30
01200	20	LMM=L2
01300	30	        DO 1 LM=1,LMM
01400		ARG=(LM-1)*SCL
01500		C=COS(ARG)
01600		S=SIN(ARG)
01700		        DO 1 LI=LIX,N,LIX
01800		J1=LI-LIX+LM
01900		J2=J1+LMX
02000		T1=X(J1)-X(J2)
02100		T2=Y(J1)-Y(J2)
02200		X(J1)=X(J1)+X(J2)
02300		Y(J1)=Y(J1)+Y(J2)
02400		X(J2)=C*T1+S*T2
02500	1	Y(J2)=C*T2-S*T1
02600		RETURN
02700		END
     

00100		SUBROUTINE RBITS(X,Y,M)
00200		DIMENSION X(512),Y(512),L(9)
00300		EQUIVALENCE(L9,L(1)),(L8,L(2)),(L7,L(3)),(L6,L(4)),
00400	     1 (L5,L(5)),(L4,L(6)),(L3,L(7)),(L2,L(8)),(L1,L(9))
00500	C       PERFORM BIT REVERSAL
00600		    DO 70 J=1,9
00700		    L(J)=1
00800		    IF (J-M) 71,71,70
00900	71	    L(J)=2**(M+1-J)
01000	70	    CONTINUE
01100		JN=1
01200		    DO 60 J1=1,L1
01300		    DO 60 J2=J1,L2,L1
01400		    DO 60 J3=J2,L3,L2
01500		    DO 60 J4=J3,L4,L3
01600		    DO 60 J5=J4,L5,L4
01700		    DO 60 J6=J5,L6,L5
01800		    DO 60 J7=J6,L7,L6
01900		    DO 60 J8=J7,L8,L7
02000		    DO 60 JR=J8,L9,L8
02100		IF (JN-JR) 61,61,62
02200	61	R=X(JN)
02300		X(JN)=X(JR)
02400		X(JR)=R
02500		FI=Y(JN)
02600		Y(JN)=Y(JR)
02700		Y(JR)=FI
02800	62	JN=JN+1
02900	60	CONTINUE
03000		RETURN
03100		END
     

00100		SUBROUTINE USCRM(X,Y,M)
00200		DIMENSION X(512),Y(512)
00300		I=2**M
00400		LO2=I/2
00500		SA=3.141593/I
00600		X(1)=2.*(X(1)+Y(1))
00700		Y(1)=0.
00800		LLL=LO2+1
00900		X(LLL)=2.*X(LLL)
01000		Y(LLL)=-2.*Y(LLL)
01100			DO 33 K=2,LO2
01200			J=I-K+2
01300			ANG=SA*(K-1)
01400			C=COS(ANG)
01500			S=SIN(ANG)
01600			AA=X(K)+X(J)
01700			BB=Y(K)-Y(J)
01800			CC=X(K)-X(J)
01900			DD=Y(K)+Y(J)
02000			XR=C*DD-S*CC
02100			XI=S*DD+C*CC
02200			X(K)=AA+XR
02300			Y(K)=BB-XI
02400			X(J)=AA-XR
02500			Y(J)=-BB-XI
02600	33	CONTINUE
02700		RETURN
02800		END